Wang-Landau sampling: Improving accuracy 



O 
(N 

o 

Q 
o 

(N 



A. A. Caparica^' and A. G. Cunha-Ncttott 
f Instituto de Fisica, Universidade Federal de Goids. CP. 131, CEP 74001-970, Goidnia, GO, Brazil 
ff Departamento de Fisica, Instituto de Giencias Exatas, 
and National Institute of Science and Technology for Complex Systems, 
Universidade Federal de Minas Gerais, CP. 702, 30123-970 Belo Horizonte, Minas Gerais, Brazil 

In this work we investigate the behavior of the microcanonical and canonical averages of the 
two-dimensional Ising model during the Wang-Landau simulation. The simulations were carried 
out using conventional Wang-Landau sampling and the 1/t scheme. Our findings reveal that the 
microcanonical average should not be accumulated during the initial modification factors / and 
outline a criterion to define this limit. We show that updating the density of states only after 
every spin-flip trials leads to a much better precision. We present a mechanism to determine 
for the given model up to what final modification factor the simulations should be carried out. 
Altogether these small adjustments lead to an improved procedure for simulations with much more 
reliable results. We compare our results with 1/t simulations. We also present an application of the 
procedure to a self-avoiding homopolymer. 
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I. INTRODUCTION 

In recent years Wang-Landau sampling (WLS')P, 
has been applied to many systems and has become a 
well-established Monte Carlo algorithm. The heuristic 
idea of the method is based on the fact that if one per- 
forms a random walk in energy space with a probability 
proportional to the reciprocal of the density of states, a 
flat histogram is generated for the energy distribution. 
Since the density of states produces huge numbers, in- 
stead of estimating g{E), the simulation is performed for 
S{E) = \ng{E), and a histogram H{E) is accumulated 
during the simulations to control the frequency of visits 
to the energy levels. At the beginning of the simulation 
we set S{E) = for all energy levels. The random walk 
in the energy space runs through all energy levels from 
with a probability 



p{E E') = min(exp [S{E) - S{E%1), (1) 

where E and E' are the energies of the current and 
the new possible configurations. Whenever a configu- 
ration is accepted we update H{E') — >■ H{E') -\- 1 and 
S{E') -> S{E') + F, where F = ln/,(/ is the so-called 
modification factor). The initial modification factor is 
taken as / = /o = e = 2.71828.... If the trial configura- 
tion is not accepted, then the current H{E) and S{E) are 
updated again. The flatness of the histogram is checked 
after a number of Monte Carlo (MC) steps and usually 
the histogram is considered flat if H{E) > 0.8(H), for 
all energies, where (H) is an average over the energies. 
If the flatness condition is fulfilled we update the modi- 
fication factor to a finer one by setting = y/Jl and 
reset the histogram H{E) = 0. Simulations are in gen- 
eral halted when In / ~ 10~*. Having in hand the density 
of states, one can calculate the canonical average of any 
thermodynamic variable as 



{X)t = 



EE{X)E9{E)e-P^ 
EE9{E)e 



-PE 



(2) 



where {X)e is the microcanonical average accumulated 
during the simulations and {3 = 1/kBT, ks is the Boltz- 
mann constant and T is the temperature. One of the 
interesting features of the method is that it can also ac- 
cess some quantities, such as the free energy and entropy, 
which are not directly available from conventional Monte 
Carlo simulations. 

As described above, the convergence of the method 
depends on both the flatness criterion and the final / 
when the simulation is interrupted, but the best choice 
of each is not obvious for each model to be studied. 

Recently some authors have asserted that although 
achieving a flat histogram is the initial motivation of the 
WLS, the flatness is not a necessary criterion to reach 
convergence . They argue that in conventional WLS 
the error saturates to a constant, while if In / decreases 
as 1/t, where i is a normalized Monte Carlo time, the er- 
ror would decrease monotonically as well. The 1/t algo- 
rithm is divided into two steps, initially the conventional 
WLS is followed, starting from S{E) ~ and then con- 
structing S{E) using a histogram updated in every new 
accepted configuration. S{E) is updated as in the con- 
ventional WLS, S{E) S{E) -I- Fi, with the initial value 
Fq = 1. After a number of moves (e.g. 1000 MC sweeps), 
we check H{E) to verify whether all the levels were vis- 
ited by the random walker at least once and then update 
Fi = Fi/2 and reset H{E) = 0. (The fiatness criterion 
is not required, even in this first stage.) Simulation is 
performed while Fi > 1/t = N/j, where j is the number 
of trial moves and N is the number of energy levels. In 
the remainder of the simulation Fi is updated every new 
configuration as Fi = 1/t up to a final chosen precision 

E final • 

The efficiency, convergence and limitations of the WLS 
has been quantitatively studied [1, [l^. In the present 
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work we perform a practical, computational study on the 
convergence and the accuracy of the method. 

In this paper we investigate the behavior of the max- 
ima of the specific heat 

C^{{E-{E)r)/T^ (3) 

and the susceptibility 

X^L'{{m^{\m\)f)/T, (4) 

where E is the energy of the configurations and m is the 
corresponding magnetization per spin during the conven- 
tional WLS and the 1/t algorithm simulations for the 
Ising model on a square lattice Q- We observe (as in 0- 
H, 01 ) that a considerable part of the conventional Wang- 
Landau simulation is not very useful because the error 
saturates. We propose some strategies to improve the 
efficiency of WLS and compare our results with exact 
calculations [ll| . Our findings lead to a new way of per- 
forming the WLS simulations. 

II. A NEW PROCEDURE FOR SIMULATIONS 
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FIG. 1. Evolution of the temperature of the extremum of the 
specific heat during the WLS, beginning from fir, for eight 
independent runs using the 80%-flatness criterion. The dots 
show where the modification factors were updated and the 
strai ght line is the result obtained using the exact data from 
Ref. [in. 



In order to test how far the simulations should go, dur- 
ing the WLS, beginning from /17, we calculate the spe- 
cific heat and the susceptibility defined in Eqs.(I3]) and 
(|4]), as well as the energy and the magnetization at some 
fixed temperatures. We use the current g{E) and from 
this time on this mean values are updated whenever we 
check the flatness of the histogram. FiglT] shows the evo- 
lution of the temperature of the maximum of the specific 
heat calculated for L = 32 for eight independent runs 
as a function of the Monte Carlo sweeps (MCS) [l3 and 
compare these results with the value obtained using the 
exact data of Ref. [ill (TdL = 32) = 2.29392979). The 
dots label the MCS when the modification factor was up- 
dated, the leftmost in each run corresponding to /17. One 
can see that around /23 = 1.1921 x 10^^ all the curves 
become stabilized in values displaced close to the exact 
value. Any further computational effort for In / < In /23 
does not lead to a better convergence. 

In order to investigate how these results are displaced 
around the exact value, we performed 100, 000 indepen- 
dent runs of WLS for L = 8 using the 80% and the 90% 
flatness criterions and built up histograms using bins of 
width 0.001. In Fig[2]we show that the histograms form 
nice Gaussians centered close, but not precisely in the 
exact value. In Fig[3|we show the same evolution for the 
temperature of the maximum of the susceptibility. One 
can see that in this case the curves do not flow to steady 
values. 

A strategy to improve the precision of the WLS is to 
update the density of states periodically (i.e., after every 
p trial configurations), instead of updating S(E) every 
spin-fiip trial. In order to investigate how this change 
affects the final result, we performed 100,000 indepen- 
dent runs {L = 8) using the 80% flatness criterion and 
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FIG. 2. Histograms of the locations of the peak of the spe- 
cific heat for the 2D Ising model during the WLS, using the 
80%- and 90%-fiatness criterions, each for 100, 000 indepen- 
dent runs, along with their best-fit Gaussians. The central 
line corresponds to the exact temperature of the maximum of 
the specific heag obtained with data from Ref. [lH . 



constructed again histograms of the locations of the peak 
of the specific heat. We tested the WLS with different 
values for p. Fig. |4| shows the Gaussian best-fits for 
p = 1 (conventional WLS), p = L and p ^ L^. The ver- 
tical line indicates the exact value using Ref. [ll|. One 
can see that the higher the values of p, the narrower the 
Gaussian curves. Defining the relative error s{X) for any 
quantity X by 
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FIG. 3. Evolution of the temperature of the extremum of the 
susceptibility during the WLS, beginning from fir, for eight 
independent runs using the 80%-flatness criterion. The dots 
show where the modification factors were updated and the 
strai ght line is the result obtained using the exact data from 
Ref. [H. 



FIG. 4. Best-fit Gaussians for the histograms of the tem- 
peratures of the peak of the specific heat for 2D Ising model 
during the WLS up to In/ = 10"*, using the 80%-flatness 
criterion, each for 100, 000 independent runs with the density 
of states being updated every p spin-flip trials. The central 
line corresponds to the exact temperature obtained with data 
from Ref. 
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we obtain the relative errors of the simulated mean values 
with respect to the result using Ref. [ll| for p = 1 , L and 
as 0.00043, 0.00019 and 0.00013, respectively. There- 
fore, we see that updating the density of states only after 
trial moves leads to more accurate results [3, [13 ■ 
In other words, adopting the Monte Carlo step (L^ trial 
moves), as it is conventional in the Metropolis algorithm, 
is also convenient in the Wang-Landau algorithm. 

In Figl5] we show the evolution of the location of the 
maximum of the heat capacity during WLS in which the 
density of states was updated only after every L^ spin- 
flip trials, beginning from /g. We see that now the curves 
flow to steady values around In / = In /13 = 1.2208 x 10 ^ 
and simulations with higher orders of the modification 
factor are unnecessary. One can see that, if on one hand 
the simulations adopting the Monte Carlo step becomes 
longer, on the other hand, the canonical averages con- 
verge to constant values much earlier. 

We now turn our attention to another important de- 
tail. What is the behavior of the microcanonical averages 
{M)e and {M'^)e during the sampling process? We 
have also evaluated the microcanonical averages during 
the simulations. In order to estimate the mean value of 
the magnetization during each flatness stage we carried 
out 1000 independent runs and calculated {M)e for each 
fi with i = 0, 1, 2, 26. In Fig|6] we show these results 
for two energy levels and see that they flow to relatively 
stable values around f-j. We therefore conclude that the 
microcanonical averages should not be accumulated be- 
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FIG. 5. Evolution of the temperature of the extremum of the 
specific heat during the WLS, beginning from /g, for eight 
independent runs. The density of states were updated after 
every L^ trial moves and the flatness criterion was 80%. The 
dots show where the modification factor was updated and the 
strai ght line is the result obtained using the exact data from 
Ref. 



fore In / = ln/7 = 7.843 x 10"^. 

In FiglT] we show the evolution of the maximum of 
the susceptibility during the simulations beginning from 
/g, updating the density of states after every L"^ spin- 
flip trials and accumulating the microcanonical averages 
only for In/ < ln/7. We observe that even for In/ = 
I11/26 ~ 10~® we do not obtain stable values like those 
of Figini However, if one takes the mean value of the 
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FIG. 6. Evolution of the microcanonical average of the mag- 
netization for 2D Ising model for L = 32 at £ = —1024 and 
— 1536 during the simulations over 1000 independent runs for 
each flatness stage. 



microcanonical averages in 24 independent runs and uses 
this result for calculating the canonical averages during 
the simulations the averages do flow to stable values, as 
shown in FiglS] This result shows that even for quantities 
that involve the magnetization the simulations can be 
carried out only up to In / = ln/13. 

The evolution of the canonical averages of the energy 
and the magnetization at a given temperature yields ev- 
idently patterns similar to those of Figl5] and FiglS] with 
the same limit modification factors. 

We would like to stress that all modification factors 
defined above using the canonical e microcanonical aver- 
ages during the simulations apply to the 2D Ising model 
and it is important to perform studies similar to those of 
Figures [SJ [5] and [5] before adopting this new procedure to 
other models to be sure on where to halt the simulations 
and where to begin accumulating the microcanonical av- 
erages. 

In view of the above observations, we propose the fol- 
lowing new procedure for simulations: 



• Instead of updating the density of states after every 
spin-fiip, we ought to update it after each Monte 
Carlo sweep; 

• WLS should be carried out only up to a In/ ~ 
^^f final defined by the canonical averages during 
the simulations; 



FIG. 7. Evolution of the temperature of the extremum of the 
susceptibility during the WLS, beginning from /g, for eight 
independent runs. The density of states was updated after 
every trial moves and the flatness criterion was 80%. The 
dots show where the modification factor was updated and the 
straight line is the result obtained using the exact data from 
Ref. [ll[ with the microcanonical average accumulated from 
ln/ = ln/7. 
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FIG. 8. Evolution of the temperature of the extremum of the 
susceptibility during the WLS, beginning from /g, for eight 
independent runs, using a common microcanonical average in 
24 independent runs. The density of states was updated after 
every L^ trial moves and the flatness criterion was 80%. The 
dots show where the modification factor was updated and the 
straight line is the result obtained using the exact data from 
Ref. [ll|] with the microcanonical average accumulated from 
ln/ = ln/7. 



III. COMPARISON WITH l/t SIMULATIONS 



• The microcanonical averages should not be accu- FigE] shows the evolution of the maxima of the specific 

mulated before In / < ln/„,cro defined by the mi- heat during the simulations using the l/t scheme, begin- 
crocanonical averages during the simulation. ning from the second stage and halting the simulations 
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when the CPU time matched up the mean time of the 
simulations of FiglS] 
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FIG. 9. Evolution of the temperature of the extremum of the 
specific heat during the 1/t simulations for eight independent 
runs beginning from the second stage. The strai ght line is the 
result obtained using the exact data from Ref Simula- 
tions were halted when the CPU time matched up the mean 
time of WLS. 

In order to compare these results we performed 
100, 000 independent runs of WLS for L = 8 up to 
In/ = ln/i3 using 80% and 90% flatness criterions 
(lFi./i3.80% and VFL./i3.90%) and built up the his- 
tograms. Next we carried out similar simulations us- 
ing the 1/t algorithm, halting the simulation when the 
CPU time matched up those of WL.fi^ {l/tgQ% and 
l/*90%)- In Fig ITU] we show the best-fit Gaussians of 
the histograms. One can see that they are not really 
centered around the exact value. The relative errors of 
the simulated mean values with respect to the result us- 
ing Ref. [ni yield 0.00041 and 0.00036, respectively, for 
WL.fi3.80% and VFL./i3.90%, and 0.0017 and 0.00081 
for 1/^80% and 1/^90%, with final Fk reaching 5.1 x 10^'' 
and 2.4 x 10""^. We see that although the widths of 
the 1/t-curves are smaller, their centers are farther apart 
from the exact value than those of WLS revealing a bi- 
ased estimation effect in the 1/t method. 



IV. FINITE-SIZE SCALING 



According to finite-size scaling theory |13l - lla | from the 
definition of the free energy one can obtain the zero field 
scaling expressions for the magnetization and the suscep- 
tibility, respectively by 



L-f/^XitL^/"). 



(6) 



(7) 
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FIG. 10. Best-fit Gaussians for the histograms of the tem- 
peratures of the peak of the specific heat for 2D Ising model 
during the WLS up to In / = 10"", using the 80%- and 90%- 
fiatness criterions, each for 100, 000 independent runs. The 
1/t simulations were carried out within the same CPU time. 
The central line corresponds to the exact temperature ob- 
tained with data from Ref. UM . 



We see that the locations of the maxima of these func- 
tions scale asymptotically as 



GqL 



(8) 



where is a quantity-dependent constant, allowing then 
the determination of T^. 

In order to compare the efficiency of the conventional 
WLS. the _//<-scheme and our procedure, we performed 
simulations with L = 32,36,40,44,48,52,56,64,72 and 
80 taking N = 24, 24, 20, 20, 20, 16, 16, 16, 12 and 12 in- 
dependent runs for each size, respectively. 

Using these scaling functions we estimated the critical 
temperature and the critical exponents /3 and 7. Taking 
a microcanonical average including all independent runs 
was important to reveal in FigI5]that for quantities that 
involve the magnetization the simulations can also be 
carried out only up to In / = In /13, but such a procedure 
does not lead to better results for the estimation of the 
canonical averages. 

Assuming = 1 we can use Eq. (jlip to determine Tc as 
the extrapolation to L — )■ 00 (L"^/" = 0) of the linear fits 
given by the locations of the maxima of the specific heat 
and the susceptibility defined by Eqs.([3])-(|3]). In FigfTTl 
we show the linear fits that converge to at L~^^'" ~ 
for conventional WLS and the new procedure, both using 
the 80% flatness criterion. The final estimate for Tc was 
taken as the mean value obtained from both fits. 

Since Tc is now estimated, we can calculate the criti- 
cal exponents /3 and 7. According to Eq.©, the maxi- 
mum of the finite- lattice susceptibility defined by Eq.Q 
is asymptotically proportional to L'*'/". In Fig[T^ we 
show these results for the conventional WLS and our 
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FIG. 11. Size dependence of the locations of the extrema in 
the specific heat and the susceptibility for conventional WLS 
(top) and using our procedure (bottom) assuming v = \. 



FIG. 12. Log-log plot of size dependence of the finite-lattice 
susceptibility at Tc{L) with 80% fiatness criterion for conven- 
tional WLS (top) and using our procedure (bottom). 



procedure, both using the 80%-criterion of flatness. In 
the vicinity of the critical temperature the magnetization 
scales as L~^/'^ . We can use Eq.® at the critical point 
to calculate the exponent /3 directly from the slope of the 
log-log graph and find /3. In FiglT3] wc show again the 
results for conventional WLS and our procedure for this 
exponent. One can sec that in all cases our procedure is 
more accurate than the conventional WLS. 

For the conventional WLS and the new procedure pro- 
posed here, simulations were carried out using 80% and 



Case 


Tc /? 7 


CPU time 


Exact 


2.2691853... 0.125 1.75 




l/t 


LIO"*^ 


2.2621(11) 0.197(14) 1.943(35) 


0.15 


5.10"^ 


2.2642(11) 0.1479(84) 1.846(18) 


0.30 


1.10"^ 


2.26848(35) 0.1297(31) 1.7833(46) 


1.51 


5.10"* 


2.26904(25) 0.1259(21) 1.7708(23) 


3.03 


1.10"* 


2.26944(11) 0.12647(94) 1.7616(17) 


15.13 


Conventional WLS 


80% 


2.26699(55) 0.1295(45) 1.7812(63) 


1.00 


90% 


2.26829(33) 0.1386(51) 1.7899(87) 


1.75 


Our procedure 


80% 


2.26934(23) 0.1270(16) 1.7631(27) 


9.78 


90% 


2.26916(12) 0.12494(68) 1.7555(32) 


22.21 



TABLE I. Finite size scaling results for the critical temper- 
ature and the critical exponents /3 and 7. The CPU times 
are expressed in terms of the time spent by the conventional 
WLS with 80%-flatness. 
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FIG. 13. Log-log plot of size dependence of the finite-lattice 
magnetization with 80% flatness criterion for conventional 
WLS at Tc = 2.26699 (top) and using our procedure at 
Tc = 2.26934 (bottom). 



90% flatness criterions and for l/t scheme the simula- 
tions were halted for In/ = lO'^, S.IO"'^, lO"'^, 5.10"^ 
and 10"^. In Table U we show the results for the l/t 
simulations, the conventional WLS and our procedure 
along with the exact values. The l/t results become ac- 
curate only when ln/^5xl0"^, and for lower values of 
In / they become worse, giving the impression that they 
are already fluctuating around the true value. The con- 
ventional WLS, displays problems of accuracy, while our 
results are adequately accurate for both 80% and 90% 
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flatness criterions. It is worthwhile mentioning that we 
have obtained high-resolution values using the 90% flat- 
ness criterion, which should be compared with the erratic 
behavior of the 1/t simulations for In / < 5 x 10~^, but 
such stringent level of flatness is difficult to apply to other 
systems |17H2(]| , resulting sometimes in non-convergence 
or even more inaccurate values. Moreover the 90% flat- 
ness criterion simulations are very time consuming. We 
conclude therefore that the widely adopted 80% flatness 
criterion is indeed the best guess, since it is applicable to 
all systems. 



V. APPLICATION TO A SELF-AVOIDING 
HOMOPOLYMER 

In this section, we apply to a homopolymer the initial 
tests to determine up to which modiflcation factor one 
should continue the WLS and from which f micro the mi- 
crocanonical averages should be accumulated. We con- 
sider a homopolymer consisting of N monomers which 
may assume any self avoiding walk (SAW) conflguration 
on a two-dimenstional lattice. 
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FIG. 14. Evolution of the microcanonical average of the 
mean square end-to-end distance for A'' = 50 at i? = —6 and 
E = —30 during the simulations over 100 independent runs 
for each flatness stage. 

Assuming that the polymer is in a bad solvent, there is 
an effective monomer- monomer attraction in addition to 
the self- avoidance constraint representing excluded vol- 
ume. For every pair of non-bonded nearest-neighbor 
monomers the energy of the polymer is reduced by e. 
The Hamiltonian for the model can be written as 



where a = 1(0) if the site i is occupied(vacant), and 
the sum is over nearest- neighbor pairs [llj. (The sum is 
understood to exclude pairs of bonded segments along 
the chain.) We used the so-called reptation or "slith- 
ering snake" move which consists of randomly adding 
a monomer from one end of the chain and removing a 
monomer from the other end, mantaining the size of the 
polymer constant. We define one Monte Carlo step as N 
attempted moves. 

In this model, besides the energy, another quantity of 
interest is the mean square end-to-end distance given by 



{R') = {[{xN-xof + {VM~Vo?]) 



T=1.5 



(10) 
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FIG. 15. Evolution of the energy at T = 1.5 during the WLS 
for eight independent runs beginning from /12. The straight 
line is the mean value obtained from 100 independent runs. 
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FIG. 16. Evolution of the mean square end-t-end distance at 
T = 4.5 during the WLS for eight independent runs beginning 
from /12. The straight Une is the mean value obtained from 
100 independent runs. 



<i-j> 



(9) As prescribed in Section II, we updated the density 

of states and the histogram only after each Monte Carlo 
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step. In order define from which modification factor to 
begin accumulating the microcanonical averages we es- 
timated {R^)e for each fi during the simulations for 
several energy levels. In FiglU] we show these results 
for E = —6 and E ~ —30 (the ground state of the 
= 50 homopolymer is Emin ~ —36). One can see 
that the microcanonical averages should be accumulated 
from f micro = fio- On the other hand, to estimate f final 
for halting the simulations, we calculated the canonical 
average of the energy and the mean square end-to-end 
distance during the simulations for five fixed tempera- 
tures, namely, T = 0.5, 1.5, 2.5, 3.5 and 4.5. In FiglTSlwe 
show the behavior of the energy at T = 1.5 and in FigfTBl 
the behavior of the mean square end-to-end distance at 
T = 4.5. All the graphs we have constructed for these 
quantities for the temperatures mentioned above have 
similar performance. We see therefore that for this model 
the simulations should be carried out up to f final = /is- 
Other more elaborate models, such as the HP model 
of protein foldin g |22| or continuous (off-lattice) models 
of polymers [23l . |24| will require a test similar to the one 
made in this section that may point to different values 
for f micro and f final- Notwithstanding that, our results 
suggest that this final modification factor will occur far 
before /26 ~ 1 + 10~^, leading to considerable CPU time 
savings. 



halting the simulations for In/ < In f final, with f final 
determined by the canonical averages during the sim- 
ulations and accumulating the microcanonical averages 
only for In / < In f micro, where f micro is found from the 
behavior of the microcanonical averages in each modifi- 
cation factor, it becomes quite accurate. Adopting the 
Monte Carlo step to update the density of states and 
delaying the start of the microcanonical averaging are 
changes that lead to a improved accuracy of the algo- 
rithm, while the proper definition of when to stop the 
simulation {f final) saves a lot of CPU time. 

It should be pointed out that a direct comparison of 
the density of states with exact calculations, although 
pictorially very impressive, is not a good test for algo- 
rithms that estimate the density of states. The canonical 
and microcanonical averages during the simulations are a 
more adequate checking parameter for convergence. An- 
other important conclusion is that no single simulation 
in particular tends to the exact value. One can obtain re- 
sults as close as possible to the exact value by increasing 
the number of independent runs. 

The great advantage of our findings is that all exist- 
ing codes using WLS can be promptly adapted to this 
new procedure just adding a few lines to the computer 
program. 



VI. CONCLUSIONS 

We have demonstrated that the conventional WLS 
presents problems of accuracy, but with very few changes 
in the implementation of the method, namely, updating 
the density of states only after each Monte Carlo step. 
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